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We study the ground state properties of spin-1 bosons in a two-dimensional optical lattice, by applying 
a variational Monte Carlo method to the S = 1 Bose-Hubbard model on a square lattice at unit filling. A 
doublon-holon binding factor introduced in the trial state provides a noticeable improvement in the vari- 
ational energy over the conventional Gutzwiller wave function and allows us to deal effectively with the 
inter-site correlations of particle densities and spins. We systematically show how spin-dependent interac- 
tions modify the superfiuid-Mott insulator transitions in the S — 1 Bose-Hubbard model due to the interplay 
between the density and spin fluctuations of bosons. Furthermore, regarding the magnetic phases in the 
Mott region, the calculated spin structure factor elucidates the emergence of nematic and ferromagnetic 
spin orders for antiferromagnetic (U2 > 0) and ferromagnetic (U2 < 0) couplings, respectively. 

KEYWORDS: Mott transition, superfluid, insulating state, nematic phase, ferromagnetism, S — 1 Bose-Hubbard 
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Recent progress in ultracold atom experiments has 
offered unprecedented opportunities for exploring fun- 
damental quantum phenomena in strongly correlated 
many-body systems that have been largely ignored in 
conventional condensed-matter physics. A prominent ex- 
ample of such phenomena is the quantum phase tran- 
sition from a superfluid (SF) to a Mott insulator (MI), 
demonstrated using cold bosons with frozen spin degrees 
of freedom trapped in optical lattices. 1 ' Furthermore, 
quantum gas microscope techniques 2-5 ) have opened the 
door to the detection and manipulation of single bosons 
at a single site level in an optical lattice, just like scan- 
ning tunneling microscopy in solid state physics. Quite 
recently, Endres et al. 6 ' used this technique to track the 
SF-MI transition in more detail, and found that corre- 
lated pairs consisting of a doubly populated site (dou- 
blon, D) and an unpopulated site (holon, H), which rep- 
resent the excitations in an MI, fundamentally determine 
the properties of the SF-MI transition, which is consis- 
tent with recent numerical studies. 7 ' 8 ' 

Now theoretical interest is naturally moving towards 
the quantum phase transitions in bosons with unfrozen 
spin degrees of freedom 9, 10 ' trapped in optical lattices. 
The doublon, which is one fragment of the elementary 
excitation in an MI, will have internal spin structures un- 
like spinless (spin-frozen) systems. Thus, we can strongly 
expect the interplay between spin correlations and the 
SF-MI transition to play key roles in the ground state of 
the systems. The simplest of such systems will be S = 1 
bosons on an optical lattice, whose properties are well 
captured by the S = 1 Bose-Hubbard model (BHM) 11 ' 
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where a JiC[ is an annihilation operator of a boson of spin 

a at the site j, hj = J2 a ^j a^j,a anQl t,Uo > 0. Here, 
a = —1,0,1, and denotes a nearest-neighbor-site 

pair; the definition of t is a half of that reported in some 
studies. The spin-dependent (last) term in eq. (1) induces 
spin mixing and enriches the physics of this model, com- 
pared with spinless models. In cold-atom systems, the 
values of Uq and U2 depend on the s-wave scattering 
wavelength characteristic of the atom species; U2 > 
(< 0) for Na (Rb) atoms. In contrast to more famil- 
iar Fermi- Hubbard models with S =1/2 where an anti- 
ferromagnetic superexchange interaction prevails in the 
strongly correlated regime, the S = 1 BHM exhibits com- 
plicated effective inter-site spin interactions that lead to 
exotic magnetic phases, 12,13 ' owing to the absence of 
the Pauli exclusion principle. Thus far, the phase dia- 
gram of this model has been studied based on mean- 
field type theories 14,15 ' including a Gutzwiller approx- 
imation (GA), 16 ' 17 ' and on density matrix renormal- 
ization group 18, 19 ' and quantum Monte Carlo (QMC) 
methods for the one-dimensional system. 20 ' 21 ' Alterna- 
tively, at the cost of the density fluctuation, an effective 
spin Hamiltonian obtained by a strong-coupling expan- 
sion 12,15,22 ' was studied to explore the magnetic struc- 
tures in the MI phase using QMC calculations in two and 
three dimensions. 23 ' 

In this letter, we provide a consistent description of the 
ground state properties and the phase transition from SF 
to MI of the S = 1 BHM on a square lattice, focusing on 
correlated pair excitations based on a variational Monte 
Carlo (VMC) approach, which is beyond conventional 
mean-field and GA techniques. We consider the simplest 
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Fig. 1. (Color online) The total energy per site of I^g) an d 
I^dh) is compared as a function of Uo/t for three system sizes. 
The VMC data of |*g> for finite L's will converge to the GA re- 
sult, the exact analytic result of I^g) f° r L — °°- The ratio U2/U0 
is fixed at 0.1, and the result is similar to that of U2 = 0. 
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Fig. 2. (Color online) (a) Condensate fraction or momentum dis- 
tribution function at k = (0, 0) and (b) density fluctuation, as a 
function of Uo/t for |* G > and of |* D H> for U 2 /U = 0.1. The inset 
in (a) shows an enlarged view near the Mott transition point. 



case of an SF-MI transition so that we restrict the parti- 
cle density p = N/N s (N: particle number, N s : the total 
number of sites) to p = 1 (unit filling) and put p = 0. 
We consider the behavior for other odd commensurate 
densities (p — 3, 5, • • • ) to be essentially the same. We 
will report even-p cases separately. 

To implement VMC calculations, we employ an occu- 
pation number representation at each site, no, ri-i). 
As a variational wave function, we use a Jastrow-type, 
|*dh) = ■Pdh-'PgI^)- Here |$) is the one-body part: 

^ = Tib ("0,1 + 4,0 + 1°)) where > 4,1 indi- 

cates the k = component of the Fourier transformation 
of cLj l . Since the total S z is well conserved in cold atom 
systems, we work in the subspace of JV Sj — 0- As an 
onsite correlation factor, the Gutzwiller projection is ex- 
tended so that it depends on the spin configurations in 
the site, 17) 

Vq = \\j(nj,i, nj >0 ,nj--i)\ni,no,ri-i)j j(m,n ,n-i\, 
3 

(2) 

where coefficients 7 are variational parameters control- 
ling rij ya , the occupation particle number of spin a on 
the site j. The dependence on spin configuration is nec- 
essary for U2 7^ 0. Since we have confirmed that the 
probability P(n) for n > 4 becomes negligible for the 
value of interest, Uo/t ^ 10, we impose a restriction, 
7(^1,^,0,^-1) = 0, on the total occupation number 
(n) at each site for rij = n^i +^0+^,-1 > 4. The D-H 
correlation factor "Pdh(?7, Tjt) used here is the same as that 
introduced in refs. 7 and 8, where r\ (rf) is a variational 
parameter (0 < 77 < 1) that controls the strength of the 
D-H binding between nearest-neighbor (lattice-diagonal) 
sites. For 77 — 1, isolated doublons and holons are prohib- 
ited. We assume 77 and 77' are independent of the onsite 
spin configuration for simplicity. It should be noted here 
that in effect a multiply populated site means a doublon 
for Uo/t ^ 20, because P(n) for n > 3 almost vanishes. 

In the VMC calculations, we first optimized the vari- 
ational parameters using a quasi-Newton method, and 
calculated the quantities for sets of model parameters 
(Uo/t, E/2/i) with several million configurations for sev- 
eral system sizes of N s = L x L sites with L = 8-24. 

We start by discussing the features of a Mott transition 
in |$dh) by comparing those obtained with a Gutzwiller 



wave function |Wq) = Vg\$) when the ratio U2/U0 is 
fixed at 0.1. Figure 1 shows how the total energy E/t is 
improved by employing |\Pdh) in the region of intermedi- 
ate correlation strength. The value of |^g) arrives at zero 
at the Brinkman-Rice transition point, 24 - 1 U$ K /t ~ 24.3. 
On the other hand, E/t of |\T/dh) is considerably less 
than that of GA, and remains finite even for large Uo/t 
values. 

In Fig. 2(a), the condensate fraction or k = (0,0) el- 
ement of the momentum distribution function, n(fe) = 
Sa(4 ofik.a), is plotted as a function of Uo/t. This 
quantity is finite in the SF phase, but should vanish as 
1/N S in the MI phase. The value of I^dh) exhibits a sud- 
den drop (for L > 16) at Uo — Uq c ~ 22i and vanishes 
as 1/V S for Uo > Uq c in the inset shown in Fig. 2(a). 
n(k = 0)/N is regarded as an order parameter of the 
Mott transition here, and so this anomaly indicates a 
first-order SF-MI transition. In fact, we confirmed that 
the small- 1 q I behavior of the density correlation func- 
tion N(q) (not shown) changes suddenly from linear to 
quadratic in momentum \q\ at Uo = Uo c - These features 
are substantially identical to those in the spinless case on 
the same D-H binding mechanism. 8 ) Next, let us look at 
the density fluctuation, ctq = (h 2 ) — (n) 2 , which is shown 
in Fig. 2(b). As is known, the erg of GA completely van- 
ishes at Uo = Uq K , which corroborates the view that 
each site is occupied by exactly one particle in the MI 
phase. On the other hand, a 2 of |\&dh) exhibits a small 
step at Uo = Uoc and remains finite for Uo > Uq c ■ This fi- 
nite density fluctuation is reflected in the small but finite 
values of -P(O) ~ f(2). In the following, we show that the 
doublon plays a crucial role for the spin structure. 

Now, we turn to the spin-dependent features caused 
by the J72-term in %. In Fig. 3(a), the optimized values 
of the D-H binding parameter 77 are plotted to recog- 
nize the effects of the £/2-term on the SF-MI transition. 
Recall that 77 controls the strength of the D-H binding 
between nearest-neighbor sites, and 7/ = 1 means that 
each doublon is tightly bound to an adjacent holon. We 
find that the SF-MI transition point Uo c /t shifts to no- 
ticeably larger values for U2/U0 = 0.3 and —0.1, whereas 
the shift for U2/U0 = 0.1 is very small. 

First, let us consider the difference in the Uo c shift be- 
tween for U2/U0 = ±0.1. Since the system is in the vicin- 
ity of the MI phase (Uo ^> t), we may restrict the Fock 
space to rij = 0,1,2 at each site; actually we confirmed 
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Fig. 3. (Color online) (a) U2/U0 dependence of optimized D-H 
binding parameter r\ as a function of Uo/t. (b),(c) show the number 
fraction of atoms with each spin-component and (d),(e) show the 
probabilities of doublons with each spin state around Uo c /t for 
U2/U0 = ±0.1. 



that rij > 3 is negligible for Uo/t ^ 20. To minimize the 
ground state energy for U 2 7^ 0, not only the number of 
doubly populated (n = 2) sites but also their spin struc- 
tures have to be optimized. To this end, it is convenient 
to employ the eigenstates of Sj, \Sj,Sj)), where Sj and 
Sj indicate the magnitude and z-component of the total 
spin at the site-j, respectively; Sj is even (odd) when 
an even (odd) number of atoms occupy the site-j. Us- 
ing ISjjSj}}, our basis set \nj,x, rijfl, for n = 2 is 
represented as 



|2,0,0> = |2,2», |1,1,0) = |2,1)), 

|0,Q»] 

f|2,0))-V2|0,0)) S 



|0,2,0) = -L(V5|2,0)) 



|1,0,1) 



1 

75 V 



(3) 



According to eq. (3), the on-site energies of n = 2 sites 
are calculated as Uq + U2 for |2, 0, 0) and |1, 1, 0), Uq for 
|0, 2, 0) , and U - U 2 for |1, 0, 1) . We also define the prob- 
ability p a = N a /N s (a = — 1, 0, 1), where N a denotes the 
total number of particles of spin a; in the present case, 
p a obeys the conditions < p a < 1 and J2 a P°< = 
Since we assume J2j $j = 0> tne relation pi = p~i holds, 
resulting in 2p\+po = 1. Using these formulae, the classi- 
cal statistical weights for the n — 2 configurations can be 
calculated as (pi) 2 for |2,0,0), 2( Pl ) 2 for 1 1,0,1), (p ) 2 
for |0, 2, 0), and 2p pi for |1, 1, 0). For U 2 ~ t, the expec- 
tation value of the C/2-term per site is similarly obtained 
as 



E 2 ( Pl )oc-U 2 pi{pi-l/2). 



(4) 



This energy has the minimum value t/2/2 at p\ = 1/4 
for U 2 < 0, and at Pl = and 1/2 for U 2 > 0. 
Actually, the VMC results in Figs. 3(b) and 3(c) show 
Pl = Ni/N a ~ 1/4 for U 2 /U = -0.1, and p x ~ 0.5 



for U 2 /Uq — 0.1. This causes imbalanced spin popula- 
tion for n = 2 sites, as shown in Figs. 3(d) and 3(e). 
Thus, we find through eq. (4) that the on-site energies 
of n = 2 sites are renormalized to Uq + U 2 /2 for U 2 < 0, 
while they are not renormalized for U 2 > 0. Since the 
onsite energies of n = 2 sites primarily govern the SF- 
MI transition, this is an appropriate explanation of the 
difference in the Uo c /t shifts between for U 2 /Uq — ±0.1. 
Finally, we point out that the degeneracy in [^-energy 
between p\ = and 0.5 for U 2 > found in eq. (4) is 
owing to the present spin-1 rotational symmetry, namely, 
\S X = 0) = (\S Z = 1) + \S Z = -l))/\/2. Therefore, we 
can generate a ground state with p\ ~ using VMC if 
we choose a certain initial condition with po 3> p\ . 

Next, we discuss the difference in the Uo c /t shifts be- 
tween for U2/U0 = 0.1 and 0.3. In the above discus- 
sion, we adopted a classical statistical weighting [E 2 
in eq. (4)], which ignores the effect of spin fluctua- 
tion caused by the t^-term, because the spin fluctu- 
ation (or singlet formation) is suppressed by the par- 
ticle hopping for small \U 2 \/Uo's. This is not the case 
for large \U2\/Uo , s. When the tVterm becomes predom- 
inant over the hopping term, the singlet state 1 0, 0) ) = 
1/V3(|0,2,0) - V2|l,0,l)) becomes si gnificant in n = 2 
sites, in order to reduce further the {/2-energy through 
the spin-exchange processes in the S^-term. In this 
case, assuming U 2 ^> t, we ignore the hopping term 
and directly diagonalize the interaction part in T-L using 
\Sj,Sj)). As a result, we find that the on-site energy of 
n = 2 sites is renormalized toward Uo — 2U2, which van- 
ishes for U2/U0 — 0.5. Thus, the value of Uo c /t is bound 
to grow rapidly as U 2 /Uq approaches 0.5, which explains 
the pronounced shift for U 2 /Uo — 0.3 in Fig. 3(a). 

Now, we study the magnetic structures and the spin 
correlations in the ground state around the SF-MI tran- 
sition. For U2/U0 > 0, the imbalance of the spin pop- 
ulations found in Fig. 3 (d) suggests that the ground 
state of the system exhibits the spin-nematic property. 
To detect this with regard to the present case, it is use- 
ful to study not only spin-correlation functions but also 
a spin-nematic parameter Q a defined as 



-E 



3^ 



(5) 



The maximum value of Q a for a single atom with S = 1 
is 1/3, which is realized in, e.g., |1, 0, 0) and |0, 0, 1), and 
also in the eigenstates of S x with S x — 0, which can 
be expressed as a coherent superposition of 1 1 , 0, 0) and 
1 0, 0, 1) . Of course, isotropic spin gives Q a = 0. Figure 
4 shows Q z for two positive values of U 2 /Uq. Reflecting 
the imbalanced spin population shown in Fig. 3 (d), large 
values of Q z are observed in the whole range of Uo/t > 0, 
that is, both in the SF and MI phases. In the MI phase, 
the single site state is not a coherent superposition of 
|1, 0, 0) and |0, 0, 1) due to the restriction of J2j S j = °- 
Thus, the spins in the MI phase exhibit a rod- like ne- 
matic structure with Q z ~ 1/3 and Q x — Q y ~ —1/6. 
Here, as mentioned in the previous paragraph, the popu- 
lation of doublon |0, 2, 0) and thus Nq increases as U 2 /Uq 
increases, which leads to lower Q z values. This is made 
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Fig. 4. (Color online) Spin nematic parameter Q z for two posi- 
tive values of U2/U0. For comparison, the GA result is also plotted, 
where, in fact, Q z cannot be determined for Uq > Uq r . 
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Fig. 5. (Color online) Example of spin structure factor for (a) 
antiferromagnetic coupling (U2 > 0) and (b) ferromagnetic cou- 
pling (U2 < 0). The phase for U /t = 10 (30) is SF (MI). The inset 
in (b) shows enlarged views of the two phases. 



manifest in that Q z for U2/U0 = 0.3 is smaller than 
that for U2/U0 = 0.1. The minima of Q z at Uo/t ~ 10 
have the same cause, i.e., in the lower Uq region, Q z 
decreases with increase in Uo/t due to an accompany- 
ing increase of U2, but, for Uo/t ^ 10, Q z enhances as 
Uq /t increases as a result of decrease of doublon density. 
In the same context, Q z is slightly reduced from the full- 
moment value 1/3 in the MI phase, where the D-H factor 
in I^dh) induces a particle density fluctuation, in con- 
trast to GA. Incidentally, the spin directions remain arbi- 
trary and indefinite for U2 = or Uq = 0, indicating that 
the point of U2/U0 = is singular. In addition, we show 
the spin structure factor, S(q) = jf- J2j,e(Sj ' S()e lq ' Ril , 
for U 2 /U = 0.1 in Fig. 5(a). The 3 constant S(q) indi- 
cates that there is no spin correlation, and the spins take 
random directions. Here, the decrease in S(q = 0) is re- 
flected in the restriction of J2j &j = 0- Thus, we find a 
nematic order is realized for U2/U0 > 0. 

For U2/U0 < 0, Q z becomes negative and its absolute 
value decreases monotonically as Uq jt increases for Uo < 
Uoc, and is almost constant ~ —1/6 in the MI phase (not 
shown) , indicating that the spins are polarized in the x-y 
plane. In Fig. 5(b), S(q) for U2/U0 = -0.1 is plotted in 
the two phases. By considering that the value of S(q = 0) 
is extremely large with almost full moment and diverges 
proportionally to N s , a ferromagnetic (FM) long-range 
order is realized in the x-y plane. The magnitude of the 
FM moment is almost constant in both the SF and MI 
phases (not shown but expected from Fig. 5(b)). The 
results reported above (for U2 s; 0) are consistent with 
those at the weak- (SF, t > U > |&2|) 17) and strong- 



(MI, t — > O) 23 ) interaction limits. 

In summary, the S — 1 Bose-Hubbard model [eq. (1)] 
on a square lattice at unit filling is studied, using a vari- 
ational Monte Carlo method. A doublon-holon binding 
factor Pdh, which is the essence of Mott transitions, not 
only improves the variational energy considerably upon 
the GA, but enables us to study the details of the spin 
structure directly even in the MI phase without resort- 
ing to an effective spin Hamiltonian. For U2 > (< 0), a 
spin nematic (ferromagnetic) phase is stabilized from SF 
to MI phases. The present results broadly support the 
phase diagrams 14-16 ^ and spin structure 17 ' 23 ) as regards 
5=1 BHM proposed in previous studies. 

Some of the numerical computations were carried out 
at the Yukawa Institute Computer Facility and at the Cy- 
berscience Center, Tohoku University. This work is sup- 
ported by a Grant-in- Aid for Scientific Research (C) and 
also by the Next Generation Supercomputing Project, 
Nanoscience Program, from MEXT of Japan. 
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